function [F_orin_data,T_axis,F_axis] = F_STFT2(t,signal,STFT_SETTINGS)
%% Short Time Fourier Transform
% author xiaoyeyimier
% date 2021/11/30
% version 2.0.1

% t:time of signal
% STFT_SETTINGS:
%  window_len     : length of window
%  STEP_T         : step of single short time fourier transform
%  Fs             : sampling frequency
%  fft_length     : length of fft
%  window_type    : 1=none 2=hamming 3=gaussian
%  gaussian_param : used only when window_type=3
len_of_t=length(t);                                                        % number of the sampling bits
Length_Of_Res=len_of_t;                                                    % remaining length after an iteration
T_number=ceil((len_of_t-STFT_SETTINGS.window_len)/STFT_SETTINGS.STEP_T)+1; % Calculating The number of sampling points on the T-axis, and the computation is compressed
index = 1;                                                                 % Calculating The times of iteration
T_axis = ((0:T_number-1)*STFT_SETTINGS.STEP_T+STFT_SETTINGS.window_len)/STFT_SETTINGS.Fs;
F_axis = (linspace(-STFT_SETTINGS.Fs/2,STFT_SETTINGS.Fs/2,STFT_SETTINGS.fft_length))/1e6;
F_orin_data = zeros(STFT_SETTINGS.fft_length,T_number);
IsDebug = false;
dt=0.1;
if(STFT_SETTINGS.window_type==1)
    win=ones(1,STFT_SETTINGS.window_len);
elseif(STFT_SETTINGS.window_type==2)
    win=my_hanning_win(STFT_SETTINGS.window_len);
else
    win=my_gaussian_win(STFT_SETTINGS.window_len,STFT_SETTINGS.gaussian_param);
end
%     QQQ=0;
timesofloop=1;
if(IsDebug)
    fig1=figure();
    set(fig1,'Position',[115.4000 51.4000 850 400]); % 设置图片位置
    aviobj = VideoWriter('C:\matlab项目\TFanalysis\Mprogramme\lib\data\test01.avi','Motion JPEG AVI');
    aviobj.Quality = 95;
    open(aviobj);
end

while Length_Of_Res >= STFT_SETTINGS.window_len
    %计算FFT一次
    signal_temp = signal(len_of_t-Length_Of_Res+1:...
        len_of_t-Length_Of_Res+STFT_SETTINGS.window_len);
    t_temp      = t(len_of_t-Length_Of_Res+1:...
        len_of_t-Length_Of_Res+STFT_SETTINGS.window_len);
    F_temp = fftshift(fft(signal_temp.*win,STFT_SETTINGS.fft_length));
    F_orin_data(:,index)=F_temp;
    index = index+1;                                                    %Number of iterations plus one
    Length_Of_Res = Length_Of_Res - STFT_SETTINGS.STEP_T;               %Update the remaining length
    data_temp = abs(F_temp)/max(abs(F_temp));
    if(IsDebug)
        peak = F_axis(data_temp==max(data_temp));
        blank_area=[peak-0.3,peak+0.3];
        data_temp=10*log10(data_temp);
        axes1 = subplot(1,2,1);
        plot(axes1,t_temp*1e3,abs(signal_temp),'k-','LineWidth',.7);
        xlabel('Time/ms','Fontname','TimesNewRoman','Fontsize',12)
        ylabel('Amptitude','Fontname','TimesNewRoman','Fontsize',12)
        axes2 = subplot(1,2,2);
        plot(axes2,F_axis,data_temp);
        hold on;
        line(axes2,[blank_area(1),blank_area(1)],[-30,5],'Color','r', 'LineWidth', 2);
        line(axes2,[blank_area(2),blank_area(2)],[-30,5],'Color','r', 'LineWidth', 2);
        line(axes2,[blank_area(1),blank_area(2)],[5,5],'Color','r', 'LineWidth', 2);
        hold off;
        axes2.YLim = [-30,5];
        xlabel('Frequency/MHz','Fontname','TimesNewRoman','Fontsize',12)
        ylabel('Normalized Amptitude/dB','Fontname','TimesNewRoman','Fontsize',12)
        frame = getframe(fig1);
        frame = frame2im(frame);
        [frame,map] = rgb2ind(frame,256);
        if(timesofloop==1)
            imwrite(frame,map,'test.gif','LoopCount',Inf,'DelayTime',dt,'ScreenSize',[800 400]);
        else
            imwrite(frame,map,'test.gif','WriteMode','append','DelayTime',dt);
        end
        F = getframe(fig1);
        writeVideo(aviobj,F);
    end
    timesofloop=timesofloop+1;
end
if IsDebug
    close(aviobj);
end
if Length_Of_Res
    zerostofill        = STFT_SETTINGS.window_len-Length_Of_Res;
    last_interval      = [signal(end-Length_Of_Res+1:end),zeros(1,zerostofill)];
    F_temp             = fftshift(fft(last_interval.*win,STFT_SETTINGS.fft_length));
    F_orin_data(:,T_number) = F_temp;                                       %Fill in the TF matrix
end
if(IsDebug),close(fig1);end
end

function GAUSSIAN_WIN = my_gaussian_win(win_len,a)
%w(t) = e^(-at^2);
if(win_len<=0)
    error('窗宽度必须为正数')
end
idx=1:win_len;
center_pos = round(win_len/2);
GAUSSIAN_WIN = 1/(2*sqrt(pi*a)).*exp(-((idx-center_pos).^2/4/a/win_len));
amp=win_len/sum(GAUSSIAN_WIN);
GAUSSIAN_WIN=GAUSSIAN_WIN*amp;
end

function HANNING_WIN = my_hanning_win(win_len)
i=1;HANNING_WIN=zeros(1,win_len);
while i<=win_len
    HANNING_WIN(i)=0.5*(1-cos(2*pi*i/(win_len-1)));
    i=i+1;
end
amp=win_len/sum(HANNING_WIN);
HANNING_WIN=HANNING_WIN*amp;
end



